Self-consistent perturbation expansion for Bose-Einstein condensates 
satisfying Goldstone's theorem and conservation laws 
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Quantum-field-theoretic descriptions of interacting condensed bosons have suffered from the lack 
of self-consistent approximation schemes satisfying Goldstone's theorem and dynamical conservation 
laws simultaneously. We present a procedure to construct such approximations systematically by 
using either an exact relation for the interaction energy or the Hugenholtz-Pines relation to express 
the thermodynamic potential in a Luttinger-Ward form. Inspection of the self-consistent perturba- 
tion expansion up to the third order with respect to the interaction shows that the two relations 
yield a unique identical result at each order, reproducing the conserving-gapless mean-field theory 
[T. Kita, J. Phys. Soc. Jpn. 74, 1891 (2005)] as the lowest-order approximation. The uniqueness 
implies that the series becomes exact when infinite terms are retained. We also derive useful ex- 
pressions for the entropy and superfluid density in terms of Green's function and a set of real-time 
dynamical equations to describe thermalization of the condensate. 
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I. INTRODUCTION 

Broken symmetry and self-consistency are among 
the most fundamental concepts in modern theoretical 
physics. The former brings a drastic change in the sys- 
tem with the appearances of the order parameter— £ and 
the corresponding Nambu-Goldstone boson^i&i The or- 
der parameter has to be determined self-consistently to- 
gether with quasiparticles responsible for the excitation. 
Thus, the latter concept is also essential for describing 
broken symmetry phases. 

Self-consistency plays crucial roles even in normal 
systems as exemplified in the Landau theory of Fermi 
liquids 5 -^ where an external perturbation produces a self- 
consistent molecular field to yield an enhanced response 
in some cases. It is also adopted commonly in various 
practical approximation schemes such as the Hartree- 
Fock theory and the density-functional theory^ with 
some infinite series incorporated in terms of the latter, 
self-consistent approximations can be far more effective 
than the simple perturbation expansion. It is worth 
pointing out that the Landau theory of Fermi liquids 
has been justified microscopically with the self-consistent 
quantum field theory^ which in turn enabled Leggett to 
extend the Landau theory to superfluid Fermi liquids. 9 

Among those self-consistent approximations is 
Baym's ^-derivable approximation with a unique prop- 
erty of obeying various dynamical conservation laws 
automatically. 10 This is certainly a character indispens- 
able for describing noncquilibrium phenomena but not 
met by the simple perturbation expansion. Indeed, the 
^-derivable approximation seems the only systematic 
microscopic scheme within the quantum field theory 
which enables us to study equilibrium and nonequilib- 
rium phenomena on an equal footing. It includes the 
Hartree-Fock theory and the Bardeen-Cooper-Schricffcr 
(BCS) theory of superconductivity as notable examples. 
Moreover, the Boltzmann equation can be derived as a 
special case of the ^-derivable approximatiom 10 ' 11 

The key functional <f> = 3>[G] above was introduced 



by Luttinger and WarcU^ as part of the exact equi- 
librium thermodynamic potential for the normal state 
in terms of the Matsubara Green's function G. Based 
on a self-consistent perturbation expansion, the expres- 
sion also provides a systematic approximation scheme 
with the desirable property of including the exact the- 
ory as a limitifii Indeed, Luttinger— subsequently used 
the Luttinger-Ward functional to obtain some general re- 
sults on normal Fermi systems such as the Fermi-surface 
sum rule anticipated by Landau^ Note that nonequilib- 
rium systems can be handled with essentially the same 
techniques as equilibrium cases by a mere change from 
the imaginary-time Matsubara contour into the real-time 
Keldysh contour ^i 1 ^ 

Thus, it would be useful to have a practical 
derivable approximation of Bose-Einstein condensates 
(BECs) with broken U(l) symmetry, where quite a few 
dynamical experiments have been carried ou t 17 i 18 ' 19 since 
the realization of the Bose-Einstein condensation with 
a trapped atomic gas in 1995.— Another key ingredi- 
ent here is the presence of a gapless excitation in the 
long wave-length limit, as first proved by Hugenholtz and 
Pinesi^i This branch, which corresponds to the Bogoli- 
ubov mode in the weak-coupling regime^ can be identi- 
fied now as the Nambu-Goldstone boson 4 of the sponta- 
neously broken U(l) symmetry. The importance of the 
two features, i.e., "conserving" and "gapless," for describ- 
ing interacting condensed bosons was already pointed 
out by Hohenberg and Martin^ in 1965 and also em- 
phasized by Griffin 2 ^ soon after the realization of the 
Bose-Einstein condensation in the trapped atomic gases. 
Despite considerable efforts, however, few systematic ap- 
proximation schemes for BEC have been known to date 
which satisfy the two fundamental properties simultane- 
ously. 

A notable exception may be the dielectric 
formalism i 25 i 26 ' 27 i 28 It is designed specifically to 
describe another important feature of interacting con- 
densed bosons that the single-particle spectrum and the 
two-particle density spectrum coincide, as first shown by 
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Gavoret and Noziere^ By incorporating local number 
conservation additionally, it has provided the gapless 
spectrum of a weakly interacting homogeneous Bose gas 
beyond the leading order. 27 However, a generalization 
of the formalism to inhomogeneous or noncquilibrium 
situations seems not straightforward when experiments 
on atomic gases are carried out with trap potentials 
dynamically 17 i 18 i 19 

Now, the main purpose of the present paper is to de- 
velop a self-consistent perturbation expansion for BEC 
satisfying Goldstone's theorem- and dynamical conser- 
vation laws simultaneously. This will be carried out by 
extending the normal-state Luttinger-Ward functional 12 
so as to obey a couple of exact relations for BEC, i.e. , that 
for the interaction energy and the Hugenholtz-Pines re- 
lation. The two relations will be shown to yield a unique 
identical result for the functional $ of BEC at least 
up to the third order in the self-consistent perturbation 
expansion. The fact indicates that the series becomes 
exact when infinite terms are retained. It also turns 
out that the expansion reproduces the conserving-gapless 
mean field theory developed earlier with a subtraction 
procedur e) 30 ' 31 as the lowest-order approximation. The 
formulation will be carried out in the coordinate space 
so that it is applicable to inhomogeneous systems such 
as those under trap potentials and with vortices. Using 
the Keldysh Green's functions, we will also extend it to 
describe nonequilibrium behaviors. Those are subjects 
with many unresolved issues^ which cannot be treated 
by other theoretical methods for BEC such as the varia- 
tional approac h 32 ! 33 and the Monte-Carlo method^ The 
whole contents here are relevant to single-particle prop- 
erties, and we are planning to discuss two-particle prop- 
erties in the near future. 

The formalism will find a wide range of applications on 
BEC. They include: (i) clarifying molecular field effects 
in condensed Bose systems corresponding to the Lan- 
dau theory of Fermi liquids; (ii) nonequilibrium phenom- 
ena of BEC such as thermalization with full account of 
the quasiparticle collisions; (iii) microscopic derivation of 
Landau two-fluid equations with definite interaction and 
temperature dependences of the viscosity coefficient, etc. 
It will also be helpful to construct a practical functional 
for BEC within the density-functional formalism, which 
still seems absent. 

It is worth pointing out finally that the Luttinger-Ward 
functional is known as the "two-particle irreducible (2PI) 
action" in the relativistic quantum-field theory ; 35 i 36 and 
the difficulties mentioned above are also encountered in 
describing its broken symmetry phases such as that of 
the 4 theory 37 ! 38 Thus, the present issue is relevant to 
a wide range of theoretical physics beyond BEC. 

This paper is organized as follows. Section [IT] sum- 
marizes exact results on an interacting Bose system in- 
cluding an expression of the thermodynamic potential for 
BEC with $. Section |nT] presents a definite procedure to 
construct <E>. It is subsequently used to obtain a first few 
series of the self-consistent perturbation expansion ex- 



plicitly. Section IIVI derives formally exact expressions of 
entropy and superfluid density in terms of Green's func- 
tion which may also be useful for their approximate eval- 
uations. Section W\ extends the formulation to describe 
nonequilibrium dynamical evolutions of BEC. Section [VTl 
summarizes the paper. We put h = = 1 throughout 
with fee the Boltzmann constant. 



II. EXACT RESULTS 

We consider identical Bose particles with mass m and 
spin described by the Hamiltonian: 



H — Hq + Hi n t , 



(1) 



with 



Hn = 



J tfWt( ri )i^( ri ), ( 2a ) 



flint =2j d3r J d?r'^{v)^{r')V{v-v')^{v')^{v). 

(2b) 

Here ip and i/j are field operators satisfying the Bose com- 
mutation relations, K% = -^V]-/! with /i as the chem- 
ical potential, and V is the interaction potential with the 
property V(r — r') = V(r' — r). Though dropped here, 
the effect of the trap potential can be included easily in 
K x . 

Let us introduce the Heisenberg representation of the 
field operators by 8 - 



V>(1) = e Tlff V(ri)e 



-TlH 



(3) 

with 1 = (n, Ti), where < t± < T 1 with T the temper- 
ature. We next express 0(1) as a sum of the condensate 
wave function \&(1) = (V'(l)) and the quasiparticle field 
0(1) as 

0(1) = *(!) + 0(1), 0(1) =#(1) + 0(1), (4) 

with (• • • ) denoting the grand-canonical average in terms 
of H. Note (0(1)) = by definition. 

Our Matsubara Green's function in the 2x2 Nambu 
space is defined in terms of and by 



G(l,2) 



TV 



(1) 
(1) 



[0(2) 0(2)]) cr 3 



G(l,2) F(l,2) 
-F(l,2) -(5(1,2) 



(5) 



where T T denotes the "time" -ordering operator— and 
(T3 is the third Pauli matrix. Every 2x2 matrix in 
the Nambu space will be distinguished with the sym- 
bol * on top of it like G. The off-diagonal elements 
F(l,2) = (T T 0(1)0(2)) and F(l,2) = (T T 0(1)0(2)) were 
introduced by Beliaev)^ 9 - which describe the pair anni- 
hilation and creation of quasiparticles inherent in BEC. 
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The factor 03 in Eq. ([5]) is usually absent in the definition 
of G*) 23 i 24 it brings an advantage that poles of G directly 
correspond to the Bogoliubov quasiparticles . 30 ' 31 

It is easily checked that the elements of G satisfy 
G*(l,2) = Giran.riTa), F(l,2) = F(2,l), G(l,2) = 
G(2,l), and F(l,2) = F*(r 2 ri, rir 2 ), with superscript 
* denoting complex conjugate. These four relations are 
expressed compactly in terms of G as 

d-3#(l,2)as = G(r2<n,ri7a), (6a) 
& 1 G*(l,2)a 1 = -G(r 1 T 2 ,r 2n ), (6b) 

where superscript t denotes Hermitian conjugate in the 
matrix algebra. Using G _1 G = 1, one can show that 
G _1 also obeys the relations of Eq. ([6]). 

Total particle number N is calculated by integrating 
(V'(l)V'(l)) over the whole space of the system, i.e., 



N= / d 3 n [¥(!)*(!) -G(l,l + )], 



(7) 



where the subscript of 1 + denotes an extra infinitesimal 
positive constant in the argument t\ to put the creation 
operator to the left for the equal-time average.— Equa- 
tion |J7J| may also be used to eliminate /x in favor of N. 
To be explicit, we will proceed by choosing (T, fx) as inde- 
pendent variables, which will be dropped in most cases. 

We summarize relevant exact results on the system be- 
low. Those of Sees. Ill K\ and iHBl have been known from 
the late 1950s. However, they often have been proved or 
reviewed only for uniform systems with different nota- 
tions. We hence provide in Appendix A detailed deriva- 
tions of those results together with that of Sec. Ill CI for 
general inhomogeneous systems so as to be compatible 
with definition (|5|) of our Green's function. 



B. Hugenholtz-Pines theorem 

As shown in Appendix IA 41 the Hugenholtz-Pines 
theorem 2 ^ can be extended to inhomogeneous systems 
as 



J d2[(V(l,2)-£(l,2)] 



" *(2) " 




' " 


.-*(2). 








(11) 



This is the condition for the excitation spectra to have a 
gapless mode as compatible with Goldstone's theorem^ 
In the homogeneous case of \1/(1) = y/n~o with no as 
the condensate density, the first row of Eq. (jlip in the 
Fourier space reduces to the familiar Hugenholtz-Pines 
relation: 21 u = E p= o — A p= o, where p = (p,ie n ) is 
the four momentum with e n = 2mrT as the Matsub- 
ara frequency (n = 0, ±1, ±2, • • • ). Note that Eq. (jTT|) is 
given here in terms of the same operator as the Dyson- 
Beliaev equation ijHJ); see the comment below Eq. (|3"9")l 
for its relevance to Goldstone's theorem. Equation JTT]) 
can also be regarded as the generalized Gross-Pitaevskii 
equatio n 40 i 41 to incorporate the quasiparticle contribu- 
tion into the self-energies. 

It is worth pointing out that our proof of Eq. (|11[) in 
Appendix IA4I has been carried out by using the gauge 
transformation relevant to the broken U(l) symmetry 
without making any specific assumptions on the struc- 
ture of the self-energies. Especially, it has removed the 
implicit supposition by Hugenholtz and Pines 21 that the 
self-energies in condensed Bose systems also be proper in 
the conventional sense. 12 



C. Interaction energy 



A. Dyson-Beliaev equation 
Equation ([5]) satisfies the Dyson-Beliaev equatiom 21 i 39 
d3[Go 1 (l,3)-E(l,3)]G(3,2) = a 8{\,2), (8) 
where Gq" 1 is defined by 



Go X (l,2) 



-& -^--& 3 K 1 ).)(}. 2). (!)) 



with (To as the 2x2 unit matrix, £ the self-energy, and 
6(1, 2) = 6( n - T 2 )6(r 1 - r 2 ). A proof of Eq. © is given 
in Appendix I A 31 



It follows from S = G 1 — G 1 and the symmetry of 



G 1 that E also obeys the relations of Eq. 
write it as 



£(1,2) 



£(1,2) A(l,2) 
-A(l,2) -E(l,2) 



Let us 



(10) 



It then follows that the elements satisfy E*(l,2) = 
Sfon.riTu), A(l,2) = A(2,l), E(l,2) = E(2,l), and 
A(l,2) = A*(r 2 Ti,nr 2 ). 



It is shown in Appendix IA 5l that the interaction energy 
(Hint), i-e., the grand-canonical average of Eq. (|2b]) . can 
be expressed in terms of the self-energies as 



(Ht, 



T 



dl / d2{2E(l,2)[*(2)*(l)-G(2,l+)] 



-A(l,2)[*(2)*(l)-F(2,l)] 
-A(l,2)[*(2)*(l)-F(2,l)]}. (12) 

This is one of the key relations indispensable below. 



D. Luttinger-Ward functional 

Luttinger and Ward 1 ^ gave an expression of the ther- 
modynamic potential Q = —TlnTre~ H / T for an inter- 
acting normal Fermi system as a functional of the self- 
energy E. Their consideration can be extended easily to 
the normal Bose system of = 0. It is more convenient 
to regard the resultant SI as a functional of G, which 
reads as 



SI = TTr[ln(-Go 1 + E) + EG] + 



(13) 
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with TrAB = JdlJ d2A(l, 2)5(2, 1+) and 

0^(1,2) = (- A _ ^5(1,2). (14) 

The quantity $ denotes contribution of all the skeleton 
diagrams in the simple perturbation expansion for SI with 
the replacement Go — > GJ^ Its functional derivative with 
respect to G yields the self-energy E as 



£(1,2) = —T~ 



S<i> 



SG(2,1) 



(15) 



It hence follows from Dyson's equation G = (Gq 1 — E) -1 
that SI is stationary with respect to G as 6Sl/6G(2,l) = 0. 

Luttinger and Ward also put Eq. (| 1 5|) into an integral 
form with the nth-order self-energy T,^ in terms of the 
interaction. To be explicit, E 1 ™) is defined as the con- 
tribution of all the nth-order skeleton diagrams in the 
simple perturbation expansion for the proper self-energy 
with the replacement Go — > G*£ Noting that there are 
2n — 1 Green's- function lines in the diagrams of £(*), Eq. 
(fl~5| can be integrated order by order into^ 



oo 1 

T V — TrE^G. 
^ 2n 



(16) 



Comparing this expression with Eq. (112p of the normal 
state (\P = 0, F = 0), we obtain a relation between the 
nth-order terms as 



= -(i/int) ( ' 

n 



(17) 



The factor 1/n is due to the extra H lrA present in the eval- 
uation of {H- m t) compared with that of $("). Hence 
the relation will hold true generally in self-consistent per- 
turbation expansions beyond the normal phase. 



E. De Dominicis-Martin theorem 

Using a series of Legendre transformations, it was 
shown by de Dominicis and Martin^ (see also Refs. 
|35| and \3(h that the thermodynamic potential SI in 
the condensed phase can be expressed as a functional 
S7[G, F, F, f ] such that 



SQ 



sn 



5G(2,1) SF(2A) 



= 0, 



= 0. 



(18) 



Thus, the exact thermodynamic potential is stationary 
with respect to variations in both the condensate wave 
function and Green's functions. Equation (|18|) general- 
izes 5il/5G(2, 1) = of the normal-state Luttinger- Ward 
functional [Eq. (p~3|) ] to condensed Bose systems. How- 
ever, no explicit SI has been known for BEC which satis- 
fies Eq. (JTBJ) via Eqs. |(5J| and (fTTj) and also includes Eq. 
(Unj) as the limit * -> 0. 



Following Eq. (| 13|) for the normal state, we now express 
S7 of the condensed phase as 



SI = -T dl d2*(l)Go 1 (l,2)*(2) 



-G^ 1 + £) + EG] + $, 



-|Tr[ln( 



(19) 



where G and G 



1 are given as Eqs. IT4]) and 
spectively, and Tr is now defined by 

£(1,2) A(l,2) 



Tr EG = dl d2 Tr 



G(2,l+ 
-F(2,l) -G(2,l_) 



-A(l,2) -£(1,2) 
F(2,l) 



(20) 



The subscript of 1_ denotes an extra infinitesimal neg- 
ative constant in n to put the creation operator to the 
left for equal-time averages, and the second Tr denotes 
the usual trace in the matrix algebra. Equation (|19p ap- 
propriately reduces to Eq. (fl3|) as W — > 0. Whereas the 
first two terms in Eq. (|19p remain finite even for the ideal 
Bose gas, 3> is made up only of contribution due to the 
interaction. This is one of the advantages for adopting 
expression ([iU]) . 

Once S7 is written as Eq. (|19[) . one can show by using 
Eq. ©, Eq. (HI]), G(l, 2) = G(2, 1), and E(l, 2) = E(2, 1) 
that condition ((18)) can be expressed equivalently with 
respect to $ as 

E(l,2) = —T " A(1,2) = 2T- 1 



9G(2,1)' 



<9F(2,1)' 
(21a) 



1 d<S> 



J d2 [£(1,2)^(2) - A(l,2)*(2)]. (21b) 

These are direct generalizations of Eq. (|15p for the normal 
state into the condensed phase. 

Now, our gapless <I>-derivable approximation denotes 
(i) constructing «3> so as to reproduce Eq. (f2"Tj) and (ii) 
determining G, and £ self-consistently with Eqs. (jSJ), 
(fTTj) . and (|21| . It will obey dynamical conservation 
law o 10 i 31 as well as Goldstone's theorem^ii thereby en- 
abling us to handle equilibrium and nonequilibrium phe- 
nomena of BEC on an equal footing with the Nambu- 
Goldstone boson. 

It is worth pointing out that expression (|19[) becomes 
exact when $ satisfies Eq. (|T71) at each order up to 
n = 00 in the self-consistent perturbation expansion. 
With £/i n t — > XHi n t in Eq. ((T|), the proof proceeds in 
exactly the same way as that of the normal stated as 
follows. First, we find with Eqs. (fTTj) and (TBJ) that 
the corresponding expression of Eq. (|19p obeys the same 
first order differential equation dQ,\/d\ = (\Hi nt }\/\ 
as the defining one S1 A = -TlnTr e -( H o+*H int )/T ^ w h ere 
(• • -)a denotes the grand-canonical average in terms of 
Ho + \H int . Second, the two expressions yield the same 
initial value Sl> = o- We hence arrive at the above conclu- 
sion. Thus, the gapless ^-derivable scheme obeying Eq. 
(TTT)) also includes the exact theory as a limit. 
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Exact relations with <f> 



A. Definite procedure for 



We now present a couple of exact relations in terms of 
$ to be satisfied in the condensed phase. Let us substi- 
tute the nth-order contribution of Eq. (|12[) into Eq. (fTT|) 
and subsequently use Eq. (|21a|) . We thereby obtain 



2n$ (n) 4 



dl / d2 



5G(^T) 



[*(2)*(1) - G(2, 1)] 



<5F(2,1) 



[*(2)*(1)- J F(2,1)] 



[*(2)*(1)- J F(2,1)] 



0. 



<5F(2,1 

Next, substitution of Eq. (flTTa) into Eq. (j21b|) yields 



(22) 



5$ 



(12 



5$ <5$ 

*(2)+2-^- -*(2) 



SG(2, 1 



<5F(2,1) 



0. (23) 



The above two equalities will play a crucial role below 
for writing $ down explicitly. 



III. CONSTRUCTING $ 

One of the basic difficulties in developing the self- 
consistent perturbation expansion for BEC may be at- 
tributed to the absence of a definite concept of "skeleton 
diagrams," which were clear in normal systems^ due 
to the appearance of finite one-particle average 'J'(l) = 
(^(1)). It brings an ambiguity as to how to count the con- 
tribution with F and W adequately in the renormalization 
process. Our approach here is to determine the contribu- 
tion to $ inherent in BEC with the exact relation of Eq. 
(|22| or Eq. (|23f so as to reproduce the Luttinger-Ward 
functional for \& — > 0, thereby avoiding the conceptual 
difficulty to define "skeleton diagrams" explicitly. The 
two relations will be shown to yield a unique result at 
each order in the self-consistent perturbation expansion 
in terms of the interaction. 

To this end, we introduce the symmetrized vertexjS 

r(°)(ll',22') = ^(r 1 -r 2 ) ( 5(T 1 -r 2 ) 

x[5(l, l')8(2,2')+ 6(1,2)8(2,1% (24) 

satisfying r(°) (IT, 22') = r(°) (22', 11') = (IT, 2'2) = 
r(°)(12', 21'). It helps us to reduce relevant Feynman di- 
agrams substantially at the expense of introducing some 
cumbersomeness in the calculation of numerical factors^ 
Our consideration below will be carried out in terms of 
topologically distinct diagrams, where G, F, and F are 
expressed in the same way as those of superconductivity - 
and r(°) is denoted by a filled circle. Following Popov, 43 
we also suppress drawing symbols for ^ and W in those 
diagrams; they can be recovered easily with the fact that 
Eq. (|2bp originally contains two pairs of creation and an- 
nihilation operators. 



Consider the nth-order contribution. The procedure 
to construct is summarized as (a)-(d) below. See 
Figs. 1 and 2 as explicit examples of relevant diagrams 
for n = 1 and 2, respectively. 

(a) Draw all the normal-state diagrams contributing to 
$("), i.e., those diagrams which appear in the Luttinger- 
Ward functional. 12 With each such diagram, associate 
the factor of the normal state. 

(b) Draw all the distinct diagrams obtained from those 
of (a) by successively changing the directions of a pair 
of incoming and outgoing arrows at each vertex. This 
enumerates all the processes where F or F is relevant 
in place of G. With each such diagram, associate an 
unknown coefficient. 

(cl) Draw all the distinct diagrams obtained from 
those of (a) and (b) by successively removing a Green's- 
function line so as to meet the condition that a further 
removal of any line from each diagram would not break it 
into two unconnected parts. The procedure incorporates 
all the processes where the condensate wave function par- 
ticipates explicitly. The latter condition guarantees that 
the self-energies obtained by Eq. (|21a[) are composed of 
connected diagrams. 

(c2) Associate an unknown coefficient with each such 
diagram, except the one consisting only of a single vertex 
in the first order, i.e., the rightmost diagram in Fig. 1, 
for which the coefficient is easily identified to be T/4. 
Indeed, the latter represents the term obtained from Eq. 
(|2b[) by replacing every field operator by its expectation 
value, i.e., the condensate wave function. 

(d) Determine the unknown coefficients of (b) and (c) 
by requiring that either Eq. (j2"2"|) or Eq. (f2"3")l be satisfied. 

It is worth pointing out that, with Eq. (|21a[) . the dia- 
grams of (cl) necessarily yields those self-energies which 
are separated into two parts by cutting a single line, 
i.e., those classified as "improper" in the conventional 
sense.— According to Eq. (|T2|) . however, we surely need 
to consider this kind of self-energy diagrams in the self- 
consistent perturbation expansion for BEC. It should also 
be mentioned that our proof of Eq. (fill , in Appendix I A 41 




00 



(a) 2 



(b)c<J> (c) C 0> 



(c)l 



FIG. 1: Diagrams contributing to Here (a)-(c) distin- 

guish three kinds of diagrams considered at different stages 
of the procedure in Sec. IIII Al and numbers and unknown 
variables Cv (v = 2b, la, lb) denote relative weights of those 
diagrams. Each coefficient should be multiplied by T/4 to 
obtain the absolute weight. 
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(a) 1 (b)c« (b)cW 

(c)c« (c)cg> (c)cg» (c)cg 

(c)c£> (c)cg> (c)cg> (c)cP) (c)cP) 

FIG. 2: Diagrams contributing to $^ 2 \ Here (a)-(c) distin- 
guish three kinds of diagrams considered at different stages of 
the procedure in Sec. IIII Al and the number 1 and unknown 

(2) 

variables cb {y = 46, • • ■ , 2e) denote relative weights of those 
diagrams. Each coefficient should be multiplied by — T/8 to 
obtain the absolute weight. 

is carried out without the implicit supposition by Hugen- 
holtz and Pines^i that the self-energies be "proper" in the 
conventional sensed Thus, there is nothing inconsistent 
on this point in our formulation. 

B. Expression of 

The first-order contribution to $ is given by the dia- 
grams of Fig. [1] They can be expressed analytically as 

= - J dlj dl' J d2 J d2'r(°>(ll',22') 

x {2G(1, l')G(2, 2') + $>F(1, 2)F(f ', 2') 
+ci 1 a ) G(l, l')*(2)*(2') + c£> [F(l, 2)*(l')* (2') 
+F(1',2')*(1)*(2)] +*(1')*(2')*(2)*(1)}, 

(25) 



where c^jj , 4a , and are unknown coefficients. 

We now require that Eq. (|22|) be satisfied, ft turns out 
that contribution of the first two diagrams in Fig. 1 van- 
ishes in the equation. This cancellation is characteristic 
of those diagrams with no condensate wave function and 
holds order by order in Eq. (|22|) . Hence relevant to Eq. 
(|22|) in Fig. I are the last three diagrams, which yield 

4 + c«=0, <£>+<#=(>, eft +2$ +2 = 0, 
respectively. We hence obtain 

$ = -1, = "4, c«=l. (26) 

Thus, the coefficients 4 have been determined uniquely. 

Alternatively, we impose Eq. (|23|) . Terms in the re- 
sultant equation can be expressed graphically by the last 
three diagrams of Fig. I with an extra symbol for the 
missing ^> at each vertex. We obtain the same equations 
as above in terms of ci], , 4a , and 41 . We hence arrive 
at Eq. ([26]) again. 

It is worth pointing out that, except for the signs of 
F and F in the definition of Eq. (0), functional (|25|) 
with coefficient ([26| is exactly identical to that of the 
conserving-gapless mean-field theor y 30 ! 31 developed ear- 
lier with a subtraction procedure. Moreover, those signs 
have been shown not to affect the physical quantities at 
all within the first order^i Thus, the mean-field theory 
has been identified here as the only self-consistent theory 
of the first order compatible with exact relations ([22l and 
(EH). 



C. Expression of $ (2) 

The second-order contribution to $ is given by the 
diagrams of Fig. [2] They are expressed analytically as 



$(2) =-? / dl '-- I rf4'r (0) (ff',22')r (0) (33 , ,44 , ){G(2,3')G(3,2 , )G(f,4')G(4,f') 
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+c%F(2', 3')F(2, 3)G(f , 4')G(4, 1') + c^F(2\ 3')F(2, 3)F(1 ', 4')F(1, 4) 

+4 2 Q } G(2, 3')G(3, 2')G(f , 4')*(4)*(l') + d$ [F{2\ 3')*(2)#(3) + F(2, 3)*(2')*(3')] G(l , 4')G(4, l') 
+cf}F{2', 3')F(2, 3)G(f , 4')*(4)*(1') + cf} [F(2', 3')¥(2)¥(3) + F(2, 3)*(2')*(3')] F(l', 4')F(1, 4) 
+4aG(2, 3')G(3, 2 / )*(l)*(4')¥(4)$(l / ) + c^ } G(2, 3')*(3)*(2')G(1, 4')*(4)*(l') 

+4? [F(2', 3')*(2)#(3) + F(2, 3)#(2')#(3')] G(f , 4')*(4)# (1') + c^F(2', 3')F(2, 3)*(l)*(4')*(4)*(l') 
+4 2 e } [^(2, 3)#(3')#(2')F(1 , 4)*(4')*(1') + F(2', 3')*(3)*(2)F(l', 4')*(4)*(1)] }. (27) 

I 

Here the first term in the curly brackets corresponds to known prefactors cl 2 "* (y = 46, • • • , 2e) are characteristic 
the normal-state process, whereas the others with un- 




• — - — • 

• ► - • 



FIG. 3: Two kinds of extra diagrams which appear in the 
calculation of Eq. !)22|) for n — 2. 

of BEC. 

We now impose Eq. on the unknown coefficients. 
It yields eleven algebraic equations originating from the 
prefactors of (i) the diagrams in the second and third 
rows of Fig. 2 and (ii) two kinds of diagrams in Fig. 3. 
They are given by 

= 4 2 J+4 = cff+cg = c£+2c2> = (28a) 

= eg? + eg? + c£> = 2$ + eg? = 2c$ + eff + 2$ 
= 2cg+<£>+4cg= Steffi, (28b) 



FIG. 5: Distinct diagrams for ' without arrows. 



The others are two kinds of diagrams in Fig. 3, for which 
there are four different ways to add a broken-line arrow. 
The corresponding equations are given by 

n o (2) i (2) r, (2) . (2) (2) . n (2) (2) . . (2) 

= 2 Aa + C 2c = 2 C2fc + 4c = C 2c + 2c 2d = C 2c + 4 4e . 

(30b) 

which replace Eq. (|28c[) . Despite the increase in the num- 
ber of equations, the solution is still given by Eq. (|29p. as 
seen easily by substituting it into Eq. (|30|) . Thus, func- 
tional (f27l) with Eq. (|29|) satisfies the Hugenholtz-Pines 
relation (|23|) besides exact relation (|22p for the interac- 
tion energy. 



n ( 2 ) i (2) , (2) (2) . (2) , (2) / 00 > 

= C 2a + C 2b + C 2c = C 2c + C 2d + 2c 2e > ( 28c ) 

respectively. Solving them, we obtain 

(2) _ (2) 
c ib ~ z ' c 4c ~ 1 > 
c (2) _ _ 4 c (2) = 2 (2) _ 4 (2) _ 2 

(2) (2) r, (2) A (2) „ (2) , , on x 

c 2a = c 2h = 2, = -4, c} 2 > = 2, = 1. (29) 

(2) 

Thus, the coefficients c v have been determined uniquely 
with Eq. (HU). 

We may alternatively require that Eq. (|23|) be satis- 
fied. Terms which appear in the calculation of Eq. (j2"3"|) 
can also be expressed graphically. They are obtained 
from the diagrams in the second and third row of Fig. 
2 and those of Fig. 3 by adding a broken-line arrow for 
the missing at a vertex in all possible ways. The in- 
sertion is topologically unique for most of the diagrams, 
e.g., those in the second row of Fig. 2; in this case the 
corresponding equation for each diagram turns out to be 
the same as that in Eq. (|28p. There are three exceptions. 
The first of them corresponds to the pair of diagrams 

(2\ 

with the coefficient c 2c in Fig. 2, which yields three dis- 
tinct diagrams of Fig. 4. Thus, 2c%] + 4c* + 2Cg& = in 
Eq. (|28b[) is now replaced by the three equations: 

= 4? + 2e£> = 2c£> + 2c g> + c£> = 4? + 4? • (30a) 
V^> 



FIG. 4: Three distinct diagrams necessary in the calculation 
of Eq. (|23p for n — 2. They are obtained from the pair of 

(2) 

diagrams with the coefficient c 2c in Fig. 2 by adding a broken- 
line arrow for the missing 'I' at a single vertex. 



D. Constructing $ ( , 

The same procedure has been used to obtain the ex- 
pression of the third-order contribution The rele- 
vant diagrams are given in Fig. [5] without arrows, and 
Fig. [6] shows additional diagrams necessary for the evalu- 
ation of Eq. (|2"2")l or Eq. (f2"3"| . The two relations have been 
checked to yield a unique result, which is summarized in 
Appendix B. 



IV. EXPRESSIONS OF ENTROPY AND 
SUPERFLUID DENSITY 

The analysis of the preceding section has clarified that 
we can generally express the thermodynamic potential of 
interacting condensed bosons as Eq. (flU)) , and $ can be 
constructed order by order uniquely so as to satisfy both 
of exact relations (f2"2"|) and (|23p . The fact implies that 
Eq. (|19p becomes exact when terms up to n — oo are 
retained in <f>. Using Eq. (|19p . we now derive formally 
exact expressions of entropy and superfluid density in 
terms of Green's function ([5]), which are also valid within 
the gapless ^-derivable approximation. We set W(l) — > 
^(ri) and 'F(l) — > ^*(ri) below as they have no explicit 
temperature dependence. 




FIG. 6: Two additional diagrams to calculate Eq. (|22[) or Eq. 
IPI for n = 3. 
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A. Entropy 



Let us expand every quantity in Eq. (fT9")) as 



G(l,2) =T 



OO 

E 

71— — OO 



(31) 



for example, with z n = 2niriT. We next transform the 
summation over z n into an integration on the complex z 
plane by using the Bosc distribution function>&i2 



/(*) 



-f(-z) 



i 



,z/T _ l ■ 



(32) 



The signs of ±/(±z) correspond to the subscripts in Eq. 

with this choice we can deform the original integra- 
tion contour encircling the imaginary axis down around 
the real axis£A2 Note -f(-z) = l + f(z). Equation fl9| 
is thereby transformed into 



f°° de » 
-V J —Trf(e){lmln[K + S(e_) - £-<t ] 

+Im£(e_)ReG(e_) + Re£(e_)ImG(e_)} + $. (33) 



Here P denotes Cauchy principal value to remove the pole 
e = of /(e) which belongs to {z n } n , Tr is defined as 
Eq. (|2"D|) without r integrals, A' = i\ (ri, r 2 ) and G(e_) = 
G(r 1 ,r 2 ;e_) with K(r 1 ,r 2 ) = a 3 KiS(ri - r 2 ) and e_ = 
e + z0_, and ReG(e_) and ImG(e_) are defined by 



ReG(e_) , gCfzI + gk), 



ImG(£- 



G(e_) - G(e + ) 
2i 



(34a) 
(34b) 



We can also express $ in terms of /(e) by using the 
Lehmann representation of G£ Just like case of the nor- 
mal system ; 44 i 45 inspection of the self-consistent pertur- 
bation series obtained in Sec. IIIII indicates that $W and 
E(") (n = 1, 2, 3) satisfy 



t 

2^ 



TrRe£(e_)ImG(e_ 



(35) 



order by order, which is apparently connected with Eq. 
(I21ap . We will proceed by assuming that Eq. holds 
generally. 

We now calculate entropy S — —(dSl/dT). Noting 
Eq. (HS|), the differentiation needs to be carried out only 
with respect to the explicit T dependence in /4i It then 
follows from Eq. (|35|) that —d$/dT exactly cancels the 
contribution of the third term in the curly bracket of Eq. 
We also use the relation df j&T = —da/de with 



a{e) = -/(e) In |/(e)| - /(-e) In |/(- e ) 



(36) 



to perform a partial integration over e. We thereby ob- 
tain 



S = P / ^a(e)Tr 



CT 



<9Re£(e_) 
de 



ImG(e_) 



-ImE(e_) 



<9ReG(e_) 
de 



(37) 



The expression is a direct extension of the normal-state 
entrop y 44 ! 45 to the system of interacting bosons with bro- 
ken U(l) symmetry. 

Adopting the mean-field approximation without the e 
dependence in the self-energy S, Eq. (|3"Tj) reduces to a 
well-known expression. To see this, let us diagonalizc 
the operator Iri = K + £ in G _1 with the Bogoliubov-de 
Gennes equation: 31 



W(ri s r a )Mr 2 ) = u v ( r i)^E u , (38) 



where E v > 0, and the eigenfunction u v (r) can be put 
into the expression: 



u u (r) v v (r) 
-<(r) -<(r) 



(39) 



with J cPr&su^, (r)o"3#„(r) = cfq8 v > u . Note that Eq. (fTT|) 
for the condensate wave function is obtained from Eq. 
(f3"5| as the limit of u v , v v — » 'J, and I£ w — * 0; hence 
there is no energy gap in the excitation energy in ac- 
cordance with Goldstone's theorem. Green's function is 
then transformed into G(ri,r 2 ;e_) = u v {vi)(s-<7n — 
Evas) -1 ^3uJ,(r 2 )<73. Substituting it into Eq. (f3"T)) and 
using the orthonormality of u v , we arrive at S = 

+ (! + WM 1 + /-)] with /* = 
The above consideration with the mean-field approxi- 
mation has exemplified that the structure of G(e_) near 
e = is directly relevant to the entropy of BEC at low 
temperatures. It also tells us that the Bogoliubov mode 
will be connected continuously to the quasiparticle mode 
which dominates the low-temperature thermal properties 
of superfluid 4 He4^ Further investigations seem required 
about the general properties of the self-energy near e = 
to elucidate this connection. 



B. Superfluid density 

We next derive an expression of the superfluid den- 
sity. Consider a homogeneous BEC where the conden- 



se) e 



lc i ri with 



sate wave function is given by \&(ri) = 
no as the condensate density. In this case, the off- 
diagonal self-energy A(l,2) acquires the spatial depen- 
dence e lq ' < - ri+r2 - ) in terms of the center-of-mass coordi- 
nate (rj + r 2 )/2. This may be realized by looking at 
the condensate contribution to A(l,2) which is given by 
V(ri - r 2 )<S(Ti - r2)*(ri)*(r 2 ). It hence follows that 
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Green's function can be expanded in terms of the basis 
function: 





-iqri 



(40) 



with V as the volume of the system, as 

G(l,2) =T^^ p ( ri )G! p (z n )^(r 2 )e- z "^-^, (41) 



with z n — 2mtiT. The momentum density (p) 
is calculated by operating — iVi to (T T ip (2)) = 
*(ri)**(r 2 ) - (3(1,2) and setting 2 = 1+ subsequently. 
Noting G(l,2) = G(2, 1) in Eq. (JSJ, the result can also 
be expressed in terms of the Nambu matrix G p (z n ) in 
Eq. ||1TJ as 



(p) = n q - 



= nq 



T 
' 2V 

T 
2V 



p+q o 

p+q 



^pTrG p (z„)i(z„). 



o-3G p (z„)i(z„) 
(42) 



Here l(z n ) is defined by 

l{Zn) = 









o *„0_ 



(43) 



and we have used n — G(l,l + ) = n with n denoting 
the particle density. We further express the term with 
G p (z„) in Eq. (gSJ as 
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^pTrG p (2„)i(z„) = ^pTri(z n )— ln[-G p 1 (z„)] 
p p 

+ ^pTri(z„)^^G p (z n ). 



It then turns out that the second term on the right-hand 
side vanishes, as can be seen by repeating the argument 
from Eq. (3.18) to Eq. (3.23) in Ref. S3 based on the 
momentum conservation of <I> at each interaction for ho- 
mogeneous systems. On the other hand, the first term 
can be transformed with the Bose distribution function / 
into an integration just below and above the real axis on 
the complex z plane. Carrying out a partial integration 
subsequently, we obtain an expression of the superfluid 
density tensor pf^ = m(d{pi) / dqj)q—o as 



( s ) X 

pi/ = mnbi 



— yv 



ds df(e) 
2ir de 







c Pi — Trlmlnl-G- 1 ^- 
dqj p 



(44) 



q=0 



This expression manifestly tells us that the structure of 
G p (e_) near e = is directly relevant to the superfluid 
density. Without e and p dependences in the self-energy, 



for example, the formula reproduces the expression given 
by Fetter— for the weakly interacting condensed Bose 
gas. It is also worth pointing out that Eq. (|44|) is identical 
in form with that of Fermi superfluids^i including the 
BCS-BEC crossover regime42 



V. REAL-TIME EQUATIONS OF MOTION 

The formulation of Sees. |TT] and IIIII can be extended 
straightforwardly to describe nonequilibrium dynamical 
evolutions of BEC. Formally, we only need to replace the 
integration contour < r < T^ 1 by the closed time-path 
contour. 15 ' 16 We here follow Keldys h 14 ' 16 to distinguish 
the forward and backward branches so that every time 
integration is limited to — oo < t < oo. The present 
transformation from equilibrium to nonequilibrium is a 
direct extension of the normal-state consideration. 45 



Let us replace 1 



V 



(rx,^) and t\ — > itl in the 



Heisenberg operators of Eq. ([3]), where superscript j = 
1, 2 distinguishes the forward (j = 1) and backward (j = 
2) branches. Using them, we next define Green's function 
in the Nambu space by 



Gy(l,2) 



Tr 



[${2?) 0(2*)]) <7 3 



0(P) 

G«(l,2) F«(l,2) 
-F y (l,2) -Gy(l,2 

with Tq as the generalized time-ordering operator 
The elements obey G*(l,2) = -G 3 -j,3-i(2, 1), 



(45) 



14,16 



and 



F«(l,2) = ^(2,1), F«(l,2) = -F 3 *-i,3-i(2,l). 
Gy (1, 2) = Gji(2, 1). The self-energy 5^/(1, 2) is defined 
similarly as Eq. (|10[) with the additional subscripts ij. 
We next introduce the 4x4 matrices (distinguished with 
" on top): 



G(l,2) 



<V(1,2) 



G n (l,2) G 12 (l,2) 
G 2 i(l,2) G 22 (l,2) 



Go 1 (1,2) 6 
6 -V(i.2) 



00 


6 " 













CT 3 = 






6 






6 





(46a) 



(46b) 



(46c) 



where Gg^l^) is given by Eq. ([9]) with d/drx —> 
—idjdt\. The self-energy matrix S is defined in the same 
way as Eq. (|46a|) . Now, the Dyson-Beliaev equation reads 
as 



/ 



d3[<V(l, 3) - (T 3 S(1, 3)a 3 ]G(3, 2) = 6(1, 2)a , (47) 
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where two 03 's originate from the path inversion for j — 
24^ Also, Eq. (fTTj) for the condensate wave function is 
replaced by 



[d2\G^(l,2)-Y / t lj (l,2)(-iy- 1 



r *(2) " 




' 0" 


L-*(2)_ 








(48) 

where we have incorporated \&(1) = ^(l 1 ) = \&(1 2 ). Note 
f (1) = here. 

Green's function in Eq. (]45p contains an extra factor i 
compared with Eq. fS"]). Taking this fact into account, the 
real-time functional $W is obtained from the equilibrium 
one of Sec. IIIII with the following modifications: (i) add 
branch indices to Green's function and the vertex [Eq. 
|)] as G — > Gij and 



rW..,(ll',22') = [-iy- 1 6 ij 6 ii ,8 jJ ,V(T 1 - T 2 )5(h i a ) 
x[«5(1,1')<5(2,2')+^(1,2')<5(2,1')],(49) 



respectively, where the factor (— l) 3 ^ 1 is due to the 
path inversion for j = 2; (ii) include summations over 
the branch indices; (iii) multiply each term by i m / 2 , 
where m is the number of condensate wave functions 
in the relevant term; (iv) multiply the resultant expres- 
sion by (— 1)™~ 1 2™/T. In steps (iii) and (iv), we have 
removed the factor T in equilibrium $ originating from 
fl = —T In Tr e~ H / T , and also considered the path change 
T\ — > it\ to reproduce the overall factor (— i) n in the nth- 
order perturbation. 

Thus, Eq. (|2"S"j) is now replaced by 

= 1 E J dl J dl 'J d2 J ^r^ y (ii',22') 

x {2GW (1, l')G^ (2, 2')+ c^F^l, 2)F i , f (l', 2') 



+4^(1, 10^(2)^(2') + ic\ x b > [F tJ (l, 2) 
x % (l')^ f (2') + F iV (1', 2')*i(l)*i(2)] 
+i 2 ^(l')^(2')*i(2)*i(l)}, (50) 

where 4^(1) = ^(l 1 ), ancl s are given as Eq. 
(f2"6")l . Real-time functionals <I>( 2 ) and <I>( 3 ) can be con- 
structed similarly from Eqs. (|27|) and (|B1[) . respectively, 
with the coefficients of Eqs. (|29|) and ()B5[> . Using 
0^.(1,2) = -G 3 - i3 -i(2,l), i^(l,2) = -^3^(2,1), 

and r<°> ,(11',22') = -^^,,3^3^, (11', 22'), one 
can show $(")* = $00. 



Accordingly, Eq. (|21a[) for the self-energies are modi- 
fied into 



£y(l,2) = (-!)*+'■ 



5$ 



(5G^(2,1) 



A ii (l,2) = -2(-l) i 



5$ 



6Fji(2,iy 



(51a) 



(51b) 



where the factor (— l) i+J is due to the two 03 's in Eq. 
(|T7|) . It follows from $(")* = $00 and the symmetries of 
G that ££(1,2) 
A y (l,2) = -A|_. !3 _, 
The functional also satisfies 



-£3-5,3-1(2,1)^(1,2) 



A^(2,l), 
(2,1), and E y (l,2) = £^(2,1). 



-^=E(- 1 ) i+i / rf2[£ ij (l,2)*,(2) 



-Ay(l, 2)^(2)], 



(52) 



which corresponds to Eq. (|21b|) . 

Equations (ITT)) , fig]), and f5"I|) form self-consistent 
equations for nonequilibrium time evolutions of BEC sat- 
isfying conservation laws and Goldstone's theorem simul- 
taneously. Approximating <I> by $ < - 1 - ) + $( 2 ) yields a non- 
vanishing collision integral, for example. Using it, we can 
describe thermalization of weakly interacting BEC mi- 
croscopically, i.e., a topic which seems not to have been 
clarified sufficiently. See, e.g., a recent book by Griffin, 
Nikuni, and Zarembai^ for the present status on this is- 
sue. 



VI. SUMMARY 

We have developed a self-consistent perturbation ex- 
pansion for BEC with broken U(l) symmetry so as to 
obey Goldstone's theorem and dynamical conservation 
laws simultaneously. First, the Luttinger-Ward thermo- 
dynamic functional for the normal stated 2 - has been ex- 
tended to a system of interacting condensed bosons as Eq. 
(fl~9|) . Next, we have presented a procedure to construct $ 
in the functional order by order with exact relations (j22|) 
and (f2"5|) in Sec. IIII Al It has been shown subsequently 
up to the third order of the self-consistent perturbation 
expansion that both of the relations yield a unique iden- 
tical result at each order as Eq. 1|25|) with Eq. (|26|) . Eq. 
([271) with Eq. J29]), and Eq. (jBTJ) with Eq. ([B5]) . This fact 
implies that the expansion converges to the exact thermo- 
dynamic potential when infinite terms are retained in $. 
Using Eq. (|19|) . we have also derived useful expressions 
for the entropy and superfluid density in terms of Green's 
function as Eqs. (|37| and (|44| . respectively. Finally, we 
have derived a set of real-time dynamical equations for 
BEC as Eqs. @7J), gHJ), and (jCT|) . 

An expansion scheme like the present one may have 
been anticipated since the work of de Dominicis and 
Martin 42 in 1964 to prove the existence of the functional 
satisfying Eq. (|18|) . However, no explicit expression for 
the functional seems to have been known to date. As 
already noted in Introduction, one of the advantages of 
the present expansion scheme over the simple perturba- 
tion expansion lies in its ability to describe nonequilib- 
rium phenomena. Another point to be mentioned is that 
it incorporates effects of the anomalous Green's function 
F(l,2) = {T T <j){\)(j){2)) more efficiently than the simple 
perturbation expansion^ This fact may be realized by 
noting that F (1,2) becomes finite with at least a single 
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interaction line in the latter scheme. Thus, nth-order 
terms with F or F in the present expansion contain ef- 
fects which show up only after the (n+l)th order in the 
simple perturbation expansion. 

We are planning to apply the present formalism to a 
wide range of systems/phenomena in BEC to elucidate 
their properties microscopically. It also remains to be 
performed to clarify two-particle correlations within the 
present formalism. 



APPENDIX A: DERIVATION OF EQS. ©, (TTTj). 
AND ([12]) 

Following the procedure sketched by Hohenberg and 
Martin^ we here derive the Dyson-Beliaev Eq. §E§ and 
the Hugenholtz-Pines relation (jTTJ) for general inhomoge- 
neous systems so as to be compatible with our definition 
[Eq. ||SJ|] of Green's function. We also prove expression 
(fT2| for the interaction energy. 



1. Time evolution operator 

Let us introduce the external perturbation ! 23 ' 35 ! 42 
ffext(ri) = f d 3 n [^( ri)77cxt (i) +^(n)»& t (l)] , 



(Al) 

where ?7 ex t(l) is periodic in t\ with the period T . The 
total Hamiltonian in this Appendix is given as a sum of 
Eqs. (HJ and J5IJ) by 



H( n ) = H + H ey[t ( n ). 



(A2) 



The extra term H ext serves as a convenient tool to derive 
various formal relations. The limit ?7 ox t — > will be taken 
once all the necessary formulas are obtained. 

We next define a time evolution operator in terms of 

= i + ]T(-l) n / dT - ■ ■ ■ 1 2 d nH{r n ) ■ ■■H(r 1 ) 

1 J Tn J Tn 



n=l 

T T exp 
T r a exp 



dnHin) 



T > Tq 



T < T 



(A3) 



where is the anti-timc-ordcring operator. Note 
U(t,t q ) -> e -( T ~ To ^ H as 77 cxt -> 0. This operator U(t, t ) 



obeys 



dU{T,T ) 

dr 



= -H(t)U(t,t ), 



(A4a) 



It also satisfies U(tq, tq) = 1 and 

U{T,Tl)U{n,T )=U(T,To). (A5) 

Equation (|A5|) is proved as follows. We see easily that 
U(t,t q ) = U(t,ti)U(ti,tq) obeys the same first-order 
differential equation with respect to r as U(t,tq). We 
also notice that the initial values at r — t\ are the same 
between the two operators, i.e., U(ti,tq) — U(ti,tq). 
We hence conclude Eq. (|A5|) . Note especially that 



U (r,ro) = U(tq,t), as can be seen easily by setting 



t = t in Eq. (|A5|) . 

2. Equations of motion 

We now introduce the Heisenberg operators: 

^(l^-ViMriMn), 

TMl^-^TOV^riMTi), 

where U(n) = U(n,0) andW-^n) = U(0,n) with < 
Ti < T . Differentiating them with respect to n and 
using Eq. ()A4|) . we obtain 

d 



(A6) 



r7 ext (l)+ [dl'V(l, 10^(10^(10^(1), (A7a) 



^-^^(1) 
%* xt (l) + / ^'7(1,10^(1)^(10^(1'). (A7b) 



with 7(1, 10 = (5(ti - r{)7(ri - ri). 

Let us define the expectation value of an arbitrary op 
erator O(l) = W-^rijo^i^^i) by 

TrT T U(/3)0(l) 



(0(1)) = 



(A8) 



with (3 = T -1 , which for 77 ext — ► reduces to the grand- 
canonical average with respect to H. We then realize 
from Eq. (|A7[) that the quantities 



*(1) = M1)>, *(1) = $(1)>, (A9) 
obey the equations of motion: 
d 



Kt *(l)=t fart (l)+» 7 (l), 



9n 



— -^ )*(!)= 77*^(1) +j?(l), 



(AlOa) 
(AlOb) 



with 



r?(l) = / dl' 7(1, 10^(10^(10^(1)), (Alia) 



dU(T,T ) 

dr 



M(t,t )H{t ). 



(A4b) 



ij(l) ee / dl' 7(1, 10(^(1)^(10^(10)- (AUb) 
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3. Dyson-Beliaev equation 

To derive Eq. ([8]), we first differentiate 9(1) in Eq. 
(|A9|) with respect to r? ext (2). Using definition (|A8|) . one 
can easily show 



-(T T *p(l)iP(2)) + 9(1)9(2) = G(1,2), 



69(1) 

^?7ext(2) 

where G(l,2) is defined by Eq. §5§ with definition (|A8 
for the expectation value. Similar calculations lead to 



69(1) 

^oxt(2) 



G(l,2), 



-^(1,2), 



<5^(1) 
tftf(l) 



% xt (2) v ' " «fy e * xt (2) 

We next introduce the self-energies by 



F(l,2), (Al2a) 
G(l,2). (A12b) 



S(l,2): 

£(1,2) 



69(2)' 

6y(i) 
59(2)' 



A(l,2) 
E(l,2) 



69(2)' 

Mm 

69(2)' 



(Al3a) 



(A13b) 



It then follows that <5?7(l)/<5?7 ex t(2), etc., can be expressed 
as 



fy(l) 
Sr] C xt(2) 



6r)(l) 69(3) 5r](l) 69(3) 



6y(i) 

6fj(l) 
6y cx t(2) 

6_m 

^xt(2) " 



69(3)6 Vcxt (2) 69(3) 6 Vcxt (2) 
d3 [E(1,3)G(3,2)-A(1,3)F(3,2)] , 

(A14a) 

d3 [-E(1,3)F(3,2) + A(1,3)G(3,2)] , 

(A14b) 

d3 [A(l, 3)G(3, 2) - E(l, 3)F(3, 2)] , 

(A14c) 

d3 [-A(l, 3)F(3, 2) + E(l, 3)G(3, 2)] . 



(A14d) 

With these preliminaries, we now differentiate Eq. (|A10j) 
with respect to r? e xt(2) or r)* Kt (2) and set 77 cxt = subse- 
quently. Using Eqs. (|A12[) and (|A14[) . one may see easily 
that the resultant four equations of motion can be writ- 
ten compactly as Eq. 



4. Hugenholtz-Pines relation 

Equation (TTTj) can be regarded as Goldstone's 
theorem^ for the broken U(l) symmetry. To derive it, 
we consider the gauge transformation: 

»fert(l)->e* x »fext(l) J ^(ri) -» e**V(ri), (A15) 



where \ 1S constant. This brings first-order changes in 
various quantities as 6t] cxt (l) = «X??oxt(l), 6r]* xt (l) = 
-iX^xt(l), 89(1) = i X 9(l), and 69(1) = -i X 9*(l). 
Collecting terms of first order in Eq. (|A10[) . we obtain 

= - 1^)69(1) - 6 Vcxt (l) - 6r)(l) 

= l x{(-^--^i)*(l)-%xt(l) 

-J d2[S(l,2)*(2)-A(l,2)f(2)]J, (Al6a) 



= Mv ^ + A V* (1) + 7?c * xt(1) 

- J d2 [A(l, 2)9(2) - E(l, 2)f (2)] J. (A16b) 

respectively, where we have performed the same trans- 
formation for 6r)(l) as Eq. (|A14ap . Noting x is arbitrary 
and comparing Eqs. (j A10[) and (|A16[) . we obtain 

r,(l) = Jd2 [E(l, 2)9(2) - A(l, 2)9(2)] , (Al7a) 

fj(l) = J d2 [-A(l, 2)9(2) + E(l, 2)9(2)] . (A17b) 

We finally substitute Eq. (|AT7|) into Eq. (|AlO|) and set 
Vcxt — 0. We thereby arrive at Eq. (fTTj) . 



5. Interaction energy 



To derive Eq. (fT2|) . we multiply Eqs. (|A7a|) and (|A7b|) 
with r] cxt = by 0(1') and 4>(1'), respectively, oper- 
ate T T , and take the thermodynamic average with Eq. 
dU and (<j)) = in mind. Noting -{T r ^-4>(1')) = 

— ■7^-(T T (f)(l)(j)(l')) + 5(1,1') and its conjugate, we can 
express the resultant equations in terms of the diagonal 
elements of Eq. §5§ as 



d2 V(l, 2)<T T ^(2)^(2)^(1)0(1 , )>, (Al8a) 



-K 1 )G(1,1') + 6(1,1') 



= / rf2U(l,2)(T T Vi(l)^(2^(2)0(l')). (A18b) 



Wc the n set 1' = 1+ and 1' = 1_ in Eqs. (|A18a|) and 
(IA18b|) . respectively. We also multiply Eqs. (|A10ap and 
(IA10b|) with next = by *(1) and *(1) from the left, 
respectively. Let us add the four equations, perform an 
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integration over 1, and multiply the result by T/4 with 
Eq. j2b}, Eq. (H, and V(l, 1') = 5fa - T[)V fa - r[) in 
mind. We thereby obtain 

d 



<9n 



^i)g(1,1')+<5(1,1') 



l'=l. 



(A19) 



We subsequently express the right-hand side of Eq. (|A19[) 
in terms of the self-energies by using Eq. ©, Eq. (jlip . 
G(l, 2) = G(2, 1), and E(l, 2) = S(2, 1). Noting that the 
subscript in G(l_,2) is effective only for T2 = t\ where 
G(l_,2) = G(l,2+), we arrive at Eq. (fT2|). 



APPENDIX B: EXPRESSION OF $ (3) 

The basic diagrams for <I>( 3 ) are given in Fig. [SJ In- 
serting arrows into them in all possible ways, we obtain 
81 distinct diagrams. The corresponding $( 3 ) may be 
expressed compactly as 



$0) = ^Tr[4GGr (0) GGr (0) GGr (0) + GGr (0) GGr (0) GGr (0) + cg^FFr (0) GGr (0) GGr (0) 

+c { 6 3 ]fGY^FGY^GGY^ + cf e FGY^FGY^GGY^ + c^FFY^FFY^GGY^ 
+c i e 3 g ) FFY^FGY^FGY^ + c^FFY^FFY^GGY^ + c { 6 f (FFY^ PGY^ FGY^ + c.c.) 
+c ( ^FFY^FFY ( -°'>FFY^ + cf} gGY^ GGY^ GGT<® + cg } gGY^ GGY^ GGY^ 

+c£\fF + c.c.)Y^GGY^GGY^ + ^(fGY^FGY^GGY^ + c . c .) + cg } (FGY^ fGY^GGY^ + c.c.) 
+4 3 / ) FfT(°)GGr(°}(3G + c.c.)r(°) + eg FGY^ FGY^ gGY^ + c^FGT^ FGY^GgY^ 
+ci 3 i ) (FgY ( -°'>FGY^GGY^ + c.c.) + cg^Gr^FGr^Gr^ + (FGY^ FgY^GGY^ + c .c.) 
+c$(fF + c.c.)Y^FFY^GGY^ + cfl{fF + c.c.)Y^ FGY^ FGY^ + ^(FFY^fGY^FGY^ + c . c .) 

+4 3 ) (FFr(°)/Fr(°)GGr( ' + c.c.) + 4 3) (/Fr(°>FGr( >FGr( > + c.c.) + 4 3) (FFr(°)/Gr(°)FGr(°) + c.c.) 

+c {3) FFY ( -°'>FFY^gGY^ + cg> (FFY^ FgY^ FGY^ + c . c .) + c£ } FFY^ FFY^ gGY^ 
+c i ^(FFY ( -°'>FGY^FgY^ + c.c.) + c^(fF + c.c.)Y^ FFY^ FFY^ + 4 3) 3Gr( 'gGr(°)GGr( ' 

+4 3) 5 Gr( )« ? Gr( )GGr(°) + 4 3) (.gGr(°)G 5 r( )GGr(°) + c.c.) + 4 3) (/>r(°)G 5 r(°)GGr( ) + c.c.) 
+4 3) (/>r(°)gGr(°)GGr(°) + c.c.) + cgc/Gr^FGr^Gr^ + c.c.) + 4 3) (/Gr( ) J FGr(°)G ff r(°) + c.c.) 
+4 3 h ) (^r(°'/Gr( )GGr(°) + c.c.) + ci-^rWi^rWGGrW + 4 3) FFr( )5Gr(°)gGr( ) 
+4t ) (FfT(°>G 5 r(°)5Gr( ) + c.c.) + 4 3) (FGr(°)F« ? r(°)gGr( ) + c.c.) + 4 3 i(i^r (0) ^GT(°>G 5 r(°) + c.c.) 
+4 3) (FGr(°)F 5 r(°) 5 Gr( ) + c.c.) + 4 3 > 5 r(°)F 5 r (0) GGr(°> + 4 3) (FFr(°)/Gr(°)/Gr(°) + c.c.) 
+4 3) (/>r(°)/Gr(°)FGr(°) + c.c.) + 4 3) (/>r(°)/Fr(°)GGr(°) + c.c.) + 4 3) (/Fr(°)FFr(°)gGr(°) + c.c.) 

+c (3 \fFY^FFY^GgY^ + c.c.) + c^ifFY^FgY^FGY^ + c.c.) + c { 4 3 v } (fFY^ FGY^ FgY^ + c.c.) 
+c (3 J,(FFY ( -°'>fFY^gGY^ + c.c.) + c^FFY^FgY^FgY^ + eg* (fFT<® fFY^ FFT<® + c.c.) 
+cfJ(FFY^FgY^FgY^ + c.c.) + 4 3 Q ) 3 ff r(°>GGr( >GGr(°> + 4 3) 5 .9r( >GGr( >GGr( } 

+4 3 7 ) (/5r (0) FGr(°)GGr(°) + c.c.) + c$(f GY^JgY^GGY^ + c.c.) + 4 3) FFr(°>GGr( >.g 5 r(°) 
+4 3) FGr( )FGr( )3«?r(°) + 4 3 >Gr(°)FGr( >« ?5 r(°) + 4 3) (//r(°)FGr(°)FGr(°) + c.c.) 
+4 3 t , ) (FFr( )//r(°)GGr( ) + c.c.) + 4 3 J(^r(°)/ ff r (0) FGr(°) + c.c.) + 4 3) (^r(°)/gr(°)FGr( ) + c.c.) 

+cfjFFY^FFY^ggY^ + cfj FFT<® FFY^ ggY^ + c (3 J gGY^ gGY^ gGY^ 
+cg ) gGY^gGY^GgY^ + (JFY^GgY^gGY^ + c.c.) + 4 3) (/Fr(°)5Gr(°)G 3 r( > + c.c.) 
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„(3) 
-3e 



(jFTWgGrWgGTW + c.c.) + c^FgT^FgT^GgT^ + c£j FgT^ FgY™ gGT<® 



3g 



+c^(fFT^fGT^FgT^ + c.c.) + (fFT^ FgT^ FgT^ + c.c.) + cf) (/FT® fFT^ fFT<® + c.c. 



-3ft, 

Here GG etc. connect adjacent two vertices with appropriately chosen arguments as 

TrGGT^GGT^GGT^ = J <&■■■[ d& G(l', 6)G(1, 6')r (0) (66', 55')G(5', 4)G(5, 4')r(°>(44', 33') 

xG(3',2)G(3,2')r(°>(22',ll'), 



(Bl) 



(B2) 



TrFFr (0) FGr (0) FGr (0) = / dl ■ ■ ■ J d6' F(l, 5)F{2, 6)r (0) (66', 55')F(6', 4')G(5', 4)r (0) (44', 33') 

xF(3', l')G(3, 2')r^ (22', 11'), (B3) 

and 5, 5, /, and / are composed of \P and 4" as 

ff (l,2) = *(l)*(2), ff(l,2) = *(l)*(2), /(1,2) = *(1)*(2), /(1,2) = *(1)*(2). (B4) 

The unknown coefficients in Eq. (jBip have been determined by Eq. ([22)) using the diagrams of Figs. [5] and [6l The 
final results are summarized as follows: 







c (3) - 


-15, 


c (3) 




- -c (3) - 6 


c (3) 

c 6f 


-12 r ( 3 )-_ r ( 3 )-S r (3) --5 


(B5a) 


c (3) 










(3) 
u 5m 


„(3) _ (8) _ 




r (3) _ fi „(3) _ (3) _ (3) _ (3) _ (3) _ fi 




c (3) 






-24, 


c 5c 




- 1 5 r (3) 

— 10, c 5 f 


= 30, 


c (3) - 12 


(B5b) 


c (3) 

L 4a 




— L 4s — 


c 4e 


= 30, 


c 4c 


_ r O) _ .(3) 


= 12, 


„(3) _ J3) _ (3) _ (3) _ ?4 (3) _ (3) _ 
L 4d ~~ L 4j ~~ °4x ~ °4t ~ L 4fc ~~ L 4y ~~ 


-15, 


c (3) 

c 46 




— L 4i — 


L 4o — 


c 4 9 — 


c 4v — 


J3) _ J3) _ 
u 4uj °47 




fi r (3) _ „0) _ (3) _ (3) _ (3) _ - 
u, c 4p — c 4z — c 4a — c 4cr — c 4fl — o, 




c 4/ 


~ c 4h 


- c (3) - 

°4m 


c 4n — 


c 4u — 




r (3) _ „(3) _ 
c 4e — L 4A — 


c (3) _ 


-6 c (3) - c (3) - c (3) - 3 


(B5c) 








■10, 


L 3h — 


C 3/ - 


r (3) _ „(3) _ 


L 3e 


- -30 r (3) - r (3) - r (3) - 1 5 r (3) - 5 


(B5d) 



It has been confirmed that Eq. (jBljl with Eq. (JB5J) also 
satisfies Eq. (ggfl . 
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